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Abstract: In this note we study the dynamics of a model recently introduced 
by one of us, that displays glassy phenomena in absence of energy barriers. Using 
an adiabatic hypothesis we derive an equation for the evolution of the energy as 
a function of time that describes extremely well the glassy behaviour observed in 
Monte Carlo simulations. 
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Despite intensive studies both on the experimental and the theoretical sides, 
the mechanisms responsible for the appearance of the glassy state are far from 
being understood. Experiments show that the relaxation time of supercooled liquids 
increases dramatically as the temperature is lowered, and the glass is formed when 
it exceeds the probing time [1]. This slowing down is related to the appearance 
of high free-energy barriers. Examples of simple models with a glassy dynamics 
resulting from high free-energy barriers have been studied in detail and can be 
found in [2, 3, 4]. Free-energy barriers can be thought as composed of an energetic 
contribution (the landscape defined by the Hamiltonian), and an entropic one related 
to the small number of paths that from a given state lead to the thermal equilibrium. 
In order to understand the importance of the entropy barriers in the formation of 
glasses, it has been recently considered a model which has not energy barriers, but 
that presents glassy phenomena due to entropy barriers [5]. 

In this note we study analytically the relaxation dynamics of the model [5] (here- 
after referred as I). 

In the model there are N distinguishable particles that can occupy N differ- 
ent states r = 1,...,N. The energy of a configuration is given by the number of 
unoccupied sites, or defining the occupation numbers of the states n r , 

N 

ff = -E^- (i) 

r=l 

We consider the sequential Metropolis dynamics such that at each Monte Carlo 
sweep a particle is chosen at random and it is moved to a random arrival state with 
probability one if the energy does not increase, and with probability exp(—/3) if the 
energy increases of a unity. It is worth noticing that the only processes decreasing 
the energy are the ones in which particles coming from departure states with = 1 
fall in arrival states with n a > 0. The energy variation of a single sweep can be zero 
or ±1. Note that in the model there are not energy barriers: the ground state can be 
always reached decreasing the energy monotonically. The ground state corresponds 
to all particles occupying one state, hence its energy per particle (or state) is — 1. 
The dynamics of this model is very slow at low temperatures because, as the energy 
decreases with time, it takes always a larger time to the system to empty one more 
state. Due to the particular rules of the dynamics it has been suggested to call it 
backgammon (BG) model and we will adopt this name hereafter 1 . 

1 We are grateful to Marc Potters for this suggestion. 
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In order to understand the relaxation dynamics of the BG model, it is useful to 
study the canonical probability distribution of the occupation numbers at equilib- 
rium. An elementary calculation, along the lines of I, shows that for large N the 
probability of the different occupation numbers factorizes, and the single number 
distribution is given by the following modified Poisson Law: 

P p (n) = (2) 

where the 'fugacity' ( is determined self-consistently by the solution of the equation: 

The internal energy is easily obtained as U — —Pp{n = 0). 

The basic observation that allows us to solve the long-time dynamics of the BG 
model is that at low temperatures the distribution (2) has a gap at low (> 0) occu- 
pation numbers. Therefore it is conceivable that in the low-temperature dynamics 
the degrees of freedom that lead to the Poisson distribution equilibrate much faster 
than the energy itself. During the dynamical evolution, the moves that change the 
energy are rare respect to the moves that lead to the equilibration of the probability 
distribution of the occupation numbers at fixed energy. Hence, at low temperatures, 
we expect the probability distribution of the occupation numbers to follow (at large 
enough times) a law of the kind eq.(2) to a high degree of approximation, with f3 
and ( related by (3) as effective time-dependent parameters. 

This observation allows us to close the dynamical equation for the evolution 
of the energy U(t), and consequently also for the occupation numbers probability 
distribution itself. It will result a theory able to predict the evolution of all the 
quantities depending just on that probability distribution. 

An adiabatic hypothesis similar to the one used here was used as an approxi- 
mation by Coolen and Sherrington to study the relaxation of the energy (and some 
other quantities) in spin glass and neural network models [6]. The present work puts 
that approximation in perspective: we believe it to be valid when only entropy, and 
not energy barriers, are the responsible for the slowing down of the dynamics. 

Let us consider a single Monte Carlo sweep, corresponding to an infinitesimal 
time interval 5t — l/N; this is specified by 3 independent random variables: 

• A state of departure d — 1, N, to be chosen with probability nd/N. 2 
2 This corresponds to the uniform choice of a particle to move. 
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• An arrival state a = 1, ...,N with a uniform probability 

• An acceptance variable x equal to 1 with probability e _/3 and zero otherwise, 
determining the acceptance of the move if it increases the energy. 

In terms of this variables the elementary energy variation can be written as 

AE(t) = -S ndil (l - S na , ) + (1 - <W)<5n ,o4,i (4) 

Observing that 1/N J2 a ^n a ,o — —U (t) we find that the average over the distributions 
of these variables reads 

[U(t + 6t)-U(t)]N = ^ = -(l + U(t))P(n = l,t)-e^U(t)(l-P(n = l,t)). (5) 

The previous expression is an exact relation. A closed equation can be obtained 
substituting the exact probability distribution P{n = l,t) at the time t with the 
modified Poisson distribution (2) for Pp*^(n = 1) with parameters (*(t) and (3*(t) 
determined self-consistently by the relations 

U(t) = -P^ t) (n = 0) ; e ,4 _ 1 = I- (6) 

The previous equation gives the expected result ^ = at equilibrium (C*(/) = 
(,[3*(t) = (3). We note that the relaxation eq.(5) is of the Markovian type; memory 
effects are not important in the relaxation of the energy. The basic mechanism of 
equilibration of the occupation number distribution on the constant energy surface 
is not expected to hold at very early stages of the dynamics (where the energy 
changes very fast). Coherently we find discrepancies with Monte Carlo simulations 
at short times, but excellent agreement for large times. The linearization of eq. (5,6) 
around the fixed point U = U eq ((3* = (3) allows to study the exponential relaxation 
time for the energy. In terms of ( as a solution of eq.(3), the relaxation time reads: 
r = [1 + (( - l)e f ][l + C - e f ]/[C 2 (l - e*)]. For large f3 it is found r ~ e^/ft 2 . 

The relaxation equation (5) is particularly simple at zero temperature. In this 
case, only processes which decrease the energy contribute to the rate of variation of 
the energy. It is easy to observe that ^ goes like exp(— I3*{t)). Because U(T*) = 
— 1 + T* at small temperatures T* this yields, 

dU ( ~ l ^ 
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Figure 1: Energy as a function of time compared to Monte Carlo data for N = 20000 
at temperatures T = 0. We do not show the results of (5,6) for very short times, 
where appreciable differences with the Monte Carlo data are observed 
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in agreement with the conjectured expression in I. The results for the decay of 
the energy at zero temperature are shown in figure 1. These have been obtained from 
the the numerical integration of equation (5) and from the Monte Carlo dynamics. 
The agreement is very good especially for large times. 

Furthermore, equation (5) is simple enough to allow for the study of the depen- 
dence of the energy on the cooling (and also heating) rate. This is shown in figure 2 
where we plot the numerical integration of equation (5) with the Monte Carlo data 
for different cooling rates (the cooling rate a is defined as the variation of temper- 
ature per unit time in the dynamical process). Similar results have been obtained 
for the dependence of the energy on the rate variation of the temperature during a 
heating process starting from the ground state configuration. All these results are 
in agreement with those presented in I and prove that hysteresis effects are very 
strong in this model. The appearance of hysteresis loops is shown in figure 2 where 
the energy is plotted as a function of the temperature during a cooling-heating pro- 
cess. Similar to what has been observed in experiments on real glasses, the enclosed 
area of the hysteresis loop decreases with the cooling-heating rates. One striking 
feature in figure 2 is that the dynamical energy in the heating process goes below 
the equilibrium one. This is a general feature observed in real glasses [7]. As far as 
we know this has never been observed in spin glasses [4]. 

We tried to interpret the decay of the energy autocorrelation function defined in 
I along the same lines, 

Cu{t > S) = -U(s)(l + U(s)) t ~ 3 (8) 

This function displays aging at low temperatures and has been used in I in or- 
der to estimate the relaxation time. Unfortunately, finding the solution for Cu{t, s) 
proves to be a much difficult task as not only the processes of emptying the states 
are relevant, but also the diffusion of 'towers' of particles from site to site con- 
tribute. Neglection of the diffusion mechanism would lead to the crude estimate 
Cu(t,s) ~ [1 + U(t)]/[1 + U(s)] which does not fit the simulation data. In order 
to disentangle the two effects we considered the micro-canonical dynamics, where 
only moves that do not variate the energy U are taken into account. In this case, 
the occupation numbers distribution thermalizes very fast while the diffusion be- 
comes uniform in time. The correlation function Cu{t, s) (which coincides with the 
canonical equilibrium one at the temperature Tjj through the equilibrium relation 
(3)) is time-translation invariant. Monte Carlo simulations in the microcanonical 
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Figure 2: Hysteresis cycles obtained integrating numerically equation (5,6) for dif- 
ferent cooling-heating rates. The cooling-heating rates are a = 3.3 x lO^MCstep^ 1 
(lines with diamonds) and a = 3.3 x 10~ 6 MCstep _1 (lines with crosses). The points 
are the Monte Carlo data for N = 20000. 
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ensemble show that Cu(t, s) is a simple exponential. We conclude that the aging in 
the function (8) is due to the interplay of the two aforementioned effects. 
Let us consider the evolution of the quantity 



M (t,s) = — 8n r (t),o8n r ( S )fi ■ 



(9) 



r 



Contributions to Mo(t,s) come from both the emptying of the states and the dif- 
fusion of the towers. In the same spirit as it has been done for the energy, we can 
write the stochastic equation for the variation of M (t, s) in an elementary move of 
a particle. We consider the zero temperature case (i.e. x — 0) and (using the same 
notation as for the energy) we get, 



where Mi(t,s) = ((^(t^i^^o) and P(n = l,t) is the average number of 
states with one particle at time t. It is clear that considering the correlations 
Mk(t,s) = (<5 nr .(t),fc^rv(s),o), h is possible to obtain in a recursive way a hierarchical 
set of first-order differential equations specifying Mj,(t, s) as a function of M^ +1 (t, s). 
Approximations can be obtained truncating the hierarchy at some order k and writ- 
ing M k+ i(t,s) ~ —P(n = k ) t)U(s). Although we did not try a systematic study, 
we expect that good approximations could be obtained truncating the hierarchy at 
a value of k large enough that P(n = k, t) is negligible for large t. Unfortunately, 
truncation at k — 1 gives very poor results for the correlation function. 

Equations (5,6) enable us to study the response function of the system defined 
as ru(t,s) = conjugated at equilibrium to the correlation function Cu(t,s) = 
{l/N 2 )J2 r ,r>(dn r (t),o8n r ,(s),o) via the fluctuation-dissipation theorem. 3 In order to 
study rjj we have emulated an aging experiment (solving equation (5)). After a 
quench from high temperature at time zero the system thermalizes during a waiting 
time t w (with (3 kept fixed), then the temperature is changed to (3 + 5(3, choosing 8f3 
small enough to be in the linear response regime. We then computed the difference 

3 Notc the difference of definition among C and Mo, in the last quantity terms of the type 
(<*n r (t),k<S<(s),o) with r r ' are absent. 



N(A t M (t, s)) = 5 nd{t)il (-5 na{t):0 5 ria{s)fi + 5 nd{s):0 ) . 
Averaging over the distributions of the variables n a and rid we obtain 



(10) 



dM (t, s) 
di 



P(n = l,t)M (t,s) + M 1 (t,s) 



(11) 
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of the energies at the time t > t w obtained with this last procedure and the one at 
P fixed 



Xu(t,t w ) 



U{t,p + 6p)-U{t,p) 
Sp 



Jt 




(12) 



In figure 3 we plot this quantity for different waiting times and temperatures as 
a function of t — t w . It is apparent that only a very weak dependence on the waiting 
time is present, i.e. Xu(t, t w ) — Xu(t—t w )- The large-time limit of Xu(t, t w ) is the one 
predicted by the fluctuation-dissipation theorem Xu(t) = —T 2 C V = dU/dp where C v 
is the equilibrium specific heat. It would be interesting to compute the correlation 
function C\j to investigate if the fluctuation-dissipation theorem is verified at all 
times. Unfortunately eqs.(5,6) do not give access to the response function conjugated 
to M (t, s) which is the correlation that exhibits aging. 

Summarizing, we have seen that some glassy features (slow relaxation and hys- 
teresis effects) of the dynamics of the BG model can be understood under the simple 
hypothesis of partial equilibrium of the relevant degrees of freedom of the system on 
the surface of constant energy. A more refined theory would be necessary to under- 
stand the aging effect in the correlation function (8) and possibly in its conjugated 
response function. 

We expect the adiabatic mechanism described in this paper to be relevant when- 
ever the entropic barriers are responsible for the slowing down of the dynamics. In 
this line, it would be interesting to study other models where this result holds and 
variants of the BG model like for instance, the Hamiltonians H = — N p-i (Z)r ^n r ,o) p 
which present a static phase transition [8]. In those similar theoretical 

analysis in the low T region to that presented in case of the BG model should be 
valid. 

Acknowledgements F.R. acknowledges Nordita for its kind hospitality dur- 
ing the realization of this work and Universidad Carlos III de Madrid for financial 
support. 
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Figure 3: Integrated response function for different temperatures T = 2.9 x 10~ 2 
(right curve), T = 4.3 x 1CT 2 (left curve) and different waiting times (ranging from 
t w = 20 to t w = 10 7 , as a function of t — t w . We see that for different temperatures 
the response function behaves roughly as a power law x(t ~ tw) — ~ a {P){t — t w ) u 
in a large range of time scales before saturating to its limiting value T 2 C V . The 
exponent v takes the value v ~ 1 nearly independent of the temperature. 
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